﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using Routines;

namespace Source
{
    public static class Mutual
    {
        private static double makeCondEntropy(long[] series, long t, double[] h1, double[] h11, double[,] h2, long partitions, long length)
        {
            // iterators
            long i, j;

            // time series values 
            // hii - normal value
            // hi  - delayed value
            long hi, hii;

            // data count 
            long count = 0;

            double hpi, hpj, pij;

            // entropy
            double cond_ent = 0.0;

            // normalization 
            double norm;

            // zero setting
            for (i = 0; i < partitions; i++)
            {
                h1[i] = h11[i] = 0;
                for (j = 0; j < partitions; j++)
                    h2[i, j] = 0;
            }

            // box count
            for (i = 0; i < length; i++)
                if (i >= t)
                {
                    hii = series[i];
                    hi = series[i - t];

                    // count points belongs to one of the box (normal time)
                    h1[hi]++;
                    // count points belongs to one of the box (deleyed time)
                    h11[hii]++;

                    // count points (normal and delayed) belongs to two boxes simultaneously
                    h2[hi, hii]++;

                    // totla count points
                    count++;
                }

            // normalizaction factor
            norm = 1.0 / (double)count;

            cond_ent = 0.0;

            // entropy evaluation
            for (i = 0; i < partitions; i++)
            {
                // probability that delayed time series value belongs to box h1[i]
                hpi = (double)(h1[i]) * norm;
                if (hpi > 0.0)
                {
                    for (j = 0; j < partitions; j++)
                    {
                        //probability that time series value belongs to box h11[j]
                        hpj = (double)(h11[j]) * norm;
                        if (hpj > 0.0)
                        {
                            // joined probability that delayed time series value belongs to box h1[i] and time series value belongs to box h11[j]
                            pij = (double)h2[i, j] * norm;
                            if (pij > 0.0)
                                cond_ent += pij * Math.Log(pij / (hpj * hpi));
                        }
                    }
                }
            }

            return cond_ent;
        }

        public static double[,] Calculate(string fileName, long partitions, long length, long exclude, long corrLength, long column, string separator)
        {
            double min, interval;

            double[] series = Routines.Routines.GetSeries(fileName, out length, exclude, column, separator);
            Routines.Routines.RescaleData(series, length, out min, out interval);

            double[] h1 = new double[partitions];
            double[] h11 = new double[partitions];
            double[,] h2 = new double[partitions, partitions];
            long[] newSeries = new long[length];

            // iterator
            long i;

            for (i = 0; i < length; i++)
                if (series[i] < 1.0)
                    // cast real time series values into discrete values from 1 to partitions
                    newSeries[i] = (long)(series[i] * (double)partitions);
                else
                    newSeries[i] = partitions - 1;

            double shannon = makeCondEntropy(newSeries, 0L, h1, h11, h2, partitions, length);
            if (corrLength >= length)
                corrLength = length - 1;


            double[,] entropy = new double[corrLength+1, 2];

            entropy[0, 0] = 0;
            entropy[0, 1] = shannon;

            long tau;

            // count entropy for all delays
            for (tau = 1; tau <= corrLength; tau++)
            {
                entropy[tau, 0] = Convert.ToDouble(tau);
                entropy[tau, 1] = makeCondEntropy(newSeries, tau, h1, h11, h2, partitions, length);
            }

            return entropy;
        }
    }
}
